

@ 4-VARIABLE MODEL: JOB GROWTH, JUNK SPREAD, SURVEY OF LOAN OFFICER AND YIELD CURVE INVERSION INDICATOR (1985:01-2011:04)@

new;
cls;

library  optmum;

load y2[109,3]=emp.txt;@different releases of job growth@
load y3[108,1]=junk.txt;@junk spread@
load y4[108,1]=indi.txt; @yield curve inversion@
load y_5[108,1]=slo_1.txt; @survey of loan officer@


y_3=zeros(108,3);
y_2=y2[1:108,.];
y_3[.,3]=y3[1:108,.];
y_3[.,2]=y4[1:108,.];

p=39;

tau=0;

t1=rows(y_2);
tau1=0;

n=t1-p+tau1-5;
c=ones(t1,1);


z_tru=y_3[1:t1,2:3]~y_5[1:t1,1]~y_2[1:t1,3];

n1=rows(z_tru);


y_for=zeros(n,26);


sig_est=zeros(n,5);
g_est=zeros(n,5);

i=1;
DO UNTIL I>N;


z1t=y_3[2+i+tau:p+i+tau,2]~y_3[2+i+tau:p+i+tau,3]~y_2[2+i+tau:p+i+tau,2];

yt=y_2[2+tau:p+i-1+tau,2]~y_3[3+tau:p+i+tau,2]~y_3[3+tau:p+i+tau,3]~y_5[3+tau:p+i+tau,1]~y_2[3+tau:p+i+tau,1];
yt_l=y_2[1+tau:p+i-2+tau,2]~y_3[2+tau:p+i-1+tau,2]~y_3[2+tau:p+i-1+tau,3]~y_5[2+tau:p+i-1+tau,1]~y_2[2+tau:p+i-1+tau,1];

z1t_l=c[i+tau:p+tau+i-2,1]~y_3[1+i+tau:p+i+tau-1,2]~y_3[1+i+tau:p+i+tau-1,3]~y_5[1+i+tau:p+i+tau-1,1]~y_2[1+i+tau:p+i+tau-1,2];


zt=y_2[1+i+tau:p+i-1+tau,2]~y_3[2+i+tau:p+i+tau,2]~y_3[2+i+tau:p+i+tau,3]~y_5[2+i+tau:p+i+tau,1]~y_2[2+i+tau:p+i+tau,2];
zt_l=y_2[i+tau:p+i-2+tau,2]~y_3[1+i+tau:p+i+tau-1,2]~y_3[1+i+tau:p+i+tau-1,3]~y_5[2+i+tau:p+i+tau,1]~y_2[1+i+tau:p+i+tau-1,2];


cov1=vcx(zt);
yy=yt';

yy=yt';
yy_l=yt_L';

yy_hat=y_2[p+i-1+tau,2]|y_3[p+i+tau,2]|y_3[p+i+tau,3]|y_5[p+i+tau,1]|y_2[p+i+tau,1];
yy_lhat=y_2[p+i-2+tau,2]|y_3[p+i-1+tau,2]|y_3[p+i-1+tau,3]|y_5[p+i-1+tau,1]|y_2[p+i-1+tau,1];

T=ROWS(YT);


prmtr_in={       -0.51517183      -0.89693832       -2.4096275       -4.3102375       -1.9324840      0.054418828       0.59570647      -0.58051305        1.6173940        1.1807367   
-0.00059321708    -0.0031136345      0.010429746        2.0425839        3.4200641     -0.027591061       0.24509054       -10.511056       0.72688704        2.8980375       0.68388057       -1.0928011      -0.26742793      0.015646561        1.0466446     -0.025927865     -0.078813716     -0.014745454   -0.00050858460 
       1.0178587     };

t=cols(yy);

@ Maximum Likelihood Estimation @


START=1;  @1951.I @



PRMTR_IN=PRMTR_IN';

                     @Initial values of the parameters @



 {xout,fout,gout,cout}=optmum(&lik_fcn,PRMTR_IN);
                     @ prmtr estimates, -log lik value, Grandient, code@

prm_fnl=trans(xout);

"Calculating Hessian..... Please be patient!!!!";
                  hessn0=hessp(&lik_fcn,xout);
                  

grdn_fnl=gradfd(&TRANS,xout);



/* Don't forget to delete the following line when you estimate the
parameters.*/





DATA = filter(xout);





F0_HAT=0|PRM_FNL[6]|PRM_FNL[7]|PRM_FNL[8]|PRM_FNL[9];
F1_HAT=(PRM_FNL[10]~PRM_FNL[11]~PRM_FNL[12]~PRM_FNL[13])|(PRM_FNL[14]~PRM_FNL[15]~PRM_FNL[16]~PRM_FNL[17])
       |(PRM_FNL[18]~PRM_FNL[19]~PRM_FNL[20]~PRM_FNL[21])|(PRM_FNL[22]~PRM_FNL[23]~PRM_FNL[24]~PRM_FNL[25]);

F_HAT=ZEROS(1,4)~1|ZEROS(4,1)~F1_HAT;
y_for[i,.]=data[t,.];


output off;




@ END OF MAIN PROGRAM @
@========================================================================@
@========================================================================@

PROC LIK_FCN(PRMTR1);

     local F,A,G,Q,R,BETA_LL,P_LL,BETA_TL,P_TL,BETA_TT,P_TT,vt,ft,
        LIK_MAT,j_iter,prmtr,sigv1,
        SIG1,SIG2,SIG3,SIG4,val,VECP_LL,F_1,F0,K,G11,G12,G13,G14,G15,
	F10,F20,F30,F40,F11,F12,F13,F14,F21,F22,F23,F24,F31,F32,F33,F34,F41,F42,F43,F44;


     PRMTR=TRANS(PRMTR1);

LOCATE 20,1;  PRMTR';


    sigv1=prmtr[1,1];

    
    sig1=prmtr[2,1];
    sig2=prmtr[3,1];
    sig3=prmtr[4,1];
    sig4=prmtr[5,1];
 
    f10=prmtr[6,1];
    f20=prmtr[7,1];
    f30=prmtr[8,1];
    f40=prmtr[9,1];

    f11=prmtr[10,1];
    f12=prmtr[11,1];
    f13=prmtr[12,1];
    f14=prmtr[13,1];  

    f21=prmtr[14,1];
    f22=prmtr[15,1];
    f23=prmtr[16,1];
    f24=prmtr[17,1];
       

    f31=prmtr[18,1];
    f32=prmtr[19,1];
    f33=prmtr[20,1];
    f34=prmtr[21,1];    

 
    f41=prmtr[22,1];
    f42=prmtr[23,1];
    f43=prmtr[24,1];
    f44=prmtr[25,1];    

    G11=prmtr[26,1];
    G12=prmtr[27,1];
    G13=prmtr[28,1];
    G14=prmtr[29,1];
    G15=prmtr[30,1];



   


@=====================================================================@


 F0=0|f10|F20|F30|f40;

 F_1= (f11~f12~f13~f14)|(f21~f22~f23~f24)|(f31~f32~f33~f34)|(f41~f42~f43~f44);

 F=(ZEROS(1,4)~1)|(ZEROS(4,1)~F_1);
	
 G=(EYE(4)~ZEROS(4,1))|(G11~G12~G13~G14~G15);

 K=(EYE(5)-G);
    
  Q=ZEROS(5,5);Q[2,2]=SIG1^2;Q[3,3]=SIG2^2;Q[4,4]=SIG3^2;Q[5,5]=SIG4^2;

    
    R=ZEROS(5,5);R[5,5]=SIGV1^2;

   
@================================================================@

 
    BETA_LL=y_2[2+tau,2]|y_3[3+tau,2]|y_3[3+tau,3]|y_5[3+tau,1]|y_2[3+tau,2];   
    
    


p_ll=eye(5)*0.000001;



     LIK_MAT=ZEROS(T,1);

J_ITER = 1;
DO UNTIL J_ITER>T;

      BETA_TL = f0+F * BETA_LL ;
      P_TL = F * P_LL * F' + Q;

      vt=yy[.,j_iter] - G * BETA_TL-K *(f)*yy_l[.,j_iter]-k*f0 ;
                  @prediction error@

      ft= G * P_TL * G' + R;       @variance of forecast error@


      BETA_TT= BETA_TL + P_TL * G' * inv(ft) * vt;
      P_TT= P_TL - P_TL * G' * inv(ft) * G * P_TL;

      LIK_MAT[J_ITER,1] =  0.5*(LN(2*pi*det(ft)) + vt'inv(ft)*vt);

      BETA_LL=BETA_TT;
      P_LL=P_TT;


J_ITER = J_ITER+1;
ENDO;

VAL = SUMC(LIK_MAT[START:T]);

LOCATE 2,20;  VAL;

RETP(VAL);
ENDP;


@>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>@

PROC FILTER(PRMTR1);

     local F,A,G,Q,R,BETA_LL,P_LL,BETA_TL,P_TL,BETA_TT,P_TT,vt,ft,y_for1,y_for2,y_for3,y_for4,y_for5,
        LIK_MAT,j_iter,prmtr,sigv1,SIG1,SIG2,SIG3,SIG4,val,VECP_LL,BETA_MAT,for_y,F_1,F0,K,
      G11,G12,G13,G14,G15,F10,F20,F30,F40,F11,F12,F13,F14,F21,F22,F23,F24,F31,F32,F33,F34,F41,F42,F43,F44,demp_filt;

     BETA_MAT=ZEROS(T,25);
      demp_filt=zeros(t,1);

     PRMTR=TRANS(PRMTR1);

LOCATE 20,1;  PRMTR';

    sigv1=prmtr[1,1];

    
    sig1=prmtr[2,1];
    sig2=prmtr[3,1];
    sig3=prmtr[4,1];
    sig4=prmtr[5,1];
 
    f10=prmtr[6,1];
    f20=prmtr[7,1];
    f30=prmtr[8,1];
    f40=prmtr[9,1];

    f11=prmtr[10,1];
    f12=prmtr[11,1];
    f13=prmtr[12,1];
    f14=prmtr[13,1];  

    f21=prmtr[14,1];
    f22=prmtr[15,1];
    f23=prmtr[16,1];
    f24=prmtr[17,1];
       

    f31=prmtr[18,1];
    f32=prmtr[19,1];
    f33=prmtr[20,1];
    f34=prmtr[21,1];    

 
    f41=prmtr[22,1];
    f42=prmtr[23,1];
    f43=prmtr[24,1];
    f44=prmtr[25,1];    

    G11=prmtr[26,1];
    G12=prmtr[27,1];
    G13=prmtr[28,1];
    G14=prmtr[29,1];
    G15=prmtr[30,1];



   


@=====================================================================@


 F0=0|f10|F20|F30|f40;

 F_1= (f11~f12~f13~f14)|(f21~f22~f23~f24)|(f31~f32~f33~f34)|(f41~f42~f43~f44);

 F=(ZEROS(1,4)~1)|(ZEROS(4,1)~F_1);
	
 G=(EYE(4)~ZEROS(4,1))|(G11~G12~G13~G14~G15);

 K=(EYE(5)-G);
    
  Q=ZEROS(5,5);Q[2,2]=SIG1^2;Q[3,3]=SIG2^2;Q[4,4]=SIG3^2;Q[5,5]=SIG4^2;

    
    R=ZEROS(5,5);R[5,5]=SIGV1^2;

   
@================================================================@



    BETA_LL=y_2[2+tau,2]|y_3[3+tau,2]|y_3[3+tau,3]|y_5[3+tau,1]|y_2[3+tau,2];   
    


p_ll=eye(5)*0.000001;

J_ITER = 1;
DO UNTIL J_ITER>T;



      BETA_TL = f0+F * BETA_LL ;
      P_TL = F * P_LL * F' + Q;

      vt=yy[.,j_iter] - G * BETA_TL-K *(f)*yy_l[.,j_iter]-k*f0 ;
                  @prediction error@

      ft= G * P_TL * G' + R;       @variance of forecast error@


      BETA_TT= BETA_TL + P_TL * G' * inv(ft) * vt;
      P_TT= P_TL - P_TL * G' * inv(ft) * G * P_TL;

   /* 
    y_for1=f0+f*beta_tt;
    y_for2=f0+f*y_for1;
    y_for3=f0+f*y_for2;
    y_for4=f0+f*y_for3;
    y_for5=f0+f*y_for4;*/

    y_for1=f0+f*beta_tt;
    y_for2=f0+f*(y_for1[1]|y_3[p+i+tau+1,2]|y_for1[3:5]);
    y_for3=f0+f*(y_for2[1]|y_3[p+i+tau+2,2]|y_for2[3:5]);
    y_for4=f0+f*(y_for3[1]|y_3[p+i+tau+3,2]|y_for3[3:5]);
    y_for5=f0+f*(y_for4[1]|y_3[p+i+tau+4,2]|y_for4[3:5]);

    demp_filt[j_iter,.]=beta_tt[5];



     BETA_MAT[J_ITER,.]=y_for1'~y_for2'~y_for3'~y_for4'~y_for5';
  

      BETA_LL=BETA_TT;
      P_LL=P_TT;

     


J_ITER = J_ITER+1;
ENDO;
RETP(BETA_MAT[START:T,.]~demp_filt);
ENDP;


@================================================================@
PROC TRANS(c0); @ constraining values of reg. coeff.@
         local c1,d1,d2;

         c1=c0;
         

         c1[1:5,.]=exp(-1*c0[1:5,.])/10;
        c1[10,.]=c0[10,1]./(1+abs(c0[10,1]));
        c1[15,.]=c0[15,1]./(1+abs(c0[15,1]));
        c1[20,.]=c0[20,1]./(1+abs(c0[20,1]));
        c1[25,.]=c0[25,1]./(1+abs(c0[25,1]));

retp(c1);
endp;

/*
Y_FOR[i,.]=beta_mat[t,.];*/
i=i+1;
endo;


tau2=2;

err_kk=z_tru[p+2+tau:n1+TAU2-6,.]-y_for[.,2:5];
err2_kk=z_tru[p+3+tau:n1+tau2-5,.]-y_for[.,7:10];
err3_kk=z_tru[p+4+tau:n1+tau2-4,.]-y_for[.,12:15];
err4_kk=z_tru[p+5+tau:n1+tau2-3,.]-y_for[.,17:20];
err5_kk=z_tru[p+6+tau:n1+tau2-2,.]-y_for[.,22:25];



mse_kk=meanc(err_kk[2:64,.].^2);

mse2_kk=meanc(err2_kk[2:64,.].^2);
mse3_kk=meanc(err3_kk[2:64,.].^2);
mse4_kk=meanc(err4_kk[2:64,.].^2);
mse5_kk=meanc(err5_kk[2:64,.].^2);






err_kk_avg=(err_kk+err2_kk+err3_kk+err4_kk+err5_kk)/5;
err1_kk_avg=(err2_kk+err3_kk+err4_kk+err5_kk)/4;

ERR_KK12_AVG=(err_kk+err2_kk)/2;
ERR_KK13_AVG=(err_kk+err2_kk+ERR3_KK)/3;

err_kk14_avg=(err_kk+err2_kk+err3_kk+err4_kk)/4;



mse_kk_avg=meanc(err_kk_avg[2:64,.].^2);
mse1_kk_avg=meanc(err1_kk_avg[2:64,.].^2);
mse_kk12_avg=meanc(err_kk12_avg[2:64,.].^2);
mse1_kk13_avg=meanc(err_kk13_avg[2:64,.].^2);
mse_kk14_avg=meanc(err_kk14_avg[2:64,.].^2);


mse_kk;
mse2_kk;
mse3_kk;
mse4_kk;
mse5_kk;


"average MSE (1-5) is ";; mse_kk_avg;
"average MSE1 (2-5)is  ";; mse1_kk_avg;
"average MSE (1-2) is ";; mse_kk12_avg;
"average MSE (1-3)is  ";; mse1_kk13_avg;
"average MSE (1-4) is ";; mse_kk14_avg;





y_for[.,5]~y_for[.,10]~y_for[.,15]~y_for[.,20]~y_for[.,25];







